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SECTION 1 
INTRODUCTION 

In colder climates the overnight frost accumulation on 
aircraft is a common occurrence. Langston [1] defined hoar frost 
as a uniform white deposit of fine crystalline texture which 
usually occurs on exposed surfaces on a cold and cloudless night 
and which is thin enough to clearly distinguish surface features 
underneath such as paint lines, markings and lettering. 
Ljungstroem [2] used the term hoar frost in the same way. In 
this paper we refer to nocturnal frost in the same sense, as the 
frost formed due to overnight blackbody cooling of surfaces. 

In some countries, including the USA, regulations do not 
permit takeoff when frost, snow or ice is adhering to transport 
aircraft. In other countries, dispatch is permitted if, in the 
judgment of the flight crew, the accumulation will not affect the 
safety of the flight. For most commercial transports, any ice 
and snow can be removed by deicing equipment. Likewise, the noc- 
turnal frost accumulated on the aircraft must also be completely 
removed prior to flight. This is usually achieved using costly 
and time-consuming glycol sprays. According to Ljungstroem [2] 
the airlines could experience a substantial savings if it were 
sufficient to clean only a small part of the wings and control 
surfaces. Brumby [3] mentioned that frost appears to have been a 
contributing factor in at least two recent takeoff accidents of 
transport aircraft. Weeks [4] "concluded that some combination 
of takeoff weight reduction and increase in takeoff speed will be 
necessary in order to maintain the normal safety margins in the 
presence of hoar frost deposits." Frost also presents a problem 
to general aviation aircraft, which have no effective means of 
frost removal. Thus a safety hazard may be present if a normal 
takeoff is attempted with frosted wings. 
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According to Brumby [3 ] , for full wing span upper surface 
roughness beginning at the leading edge and extending varying 
distances aft, the typical effects are a reduction in maximum 
lift, a reduction in the angle of attack at which stall occurs, 
and a rapid post-stall drag increase. The effects become more 
adverse as the size and chordwise extent of the roughness 
increase. They may also be accompanied by a reduction in lift at 
a given angle of attack and by an increase in the wing parasite 
drag. Controversy concerning performance degradation with frost 
or snow on the wings may have arisen from tests on military 
aircraft showing that frost appeared to cause no degradation in 
takeoff performance. The tests, according to Langston, were 
directed only at establishing if the particular aircraft would 
take off at the handbook speeds. No attempt was made to deter- 
mine how much the stall margin had been reduced by the frost. 
Langston [1] , in citing this conclusion, added the observations 
that the hoar frost layer was up to 6 or 7 mm thick in rare 
cases. Weeks [4] stated that a frost deposit on the aircraft was 
typically 0.25 mm thick and an extremely severe deposit was about 
1 mm thick. One wonders if the rare frost observed by Langston 
is actually rime frost due to overnight drizzles, which may pro- 
duce a thickness of 6 to 7 mm. Thus, Langston [1] recommended 
that the frost should continue to be cleared off until better 
knowledge was obtained about frost and its effects on takeoff 
performance. According to Langston[l], once the aircraft is air- 
borne, the frost usually dissipates rapidly due to aerodynamic 
heating and a lower atmospheric humidity. He suggested research 
areas to clear up the confusion concerning the safety of takeoff 
with a frosted wing. 

This paper addresses several research areas suggested by 
Langston. In the second and third sections of this paper, the 
nocturnal frost formation model is described for the prediction 
of frost density and thickness as a function of time and varying 
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atmospheric conditions. This model was adapted and modified from 
the basic frost formation model developed by Dietenberger [5] . 
Comparison with a set of nocturnal frost formation experiments on 
a slightly inclined plate performed by UDRI in the 1980-81 winter 
is discussed in the fourth section. In the fifth section, an 
empirical relationship is developed for the lift, drag, and angle 
of attack penalties that would occur on a frost-coated airfoil as 
a function of the chord Reynold's number, the average roughness 
height to chord ratio, and the percent of frost-coated wing 
section. The data for the curve fittings was taken from 
Ljungstroem [2], Brumby [3], and others. The approach taken here 
is also a generalization of Brumby's method for the lift 
penalties. To a certain degree the results obtained can be 
extended to other wing section types and Reynold's number. Also, 
in this section, takeoff modifications due to aircraft roughness 
are suggested. In the last section, overnight atmospheric 
profiles are forecasted from the meteorological measurements at 5 
p.m. according to equations obtained from Bagdonas, et.al. [6]. 
The nocturnal frost formation model was then used for a wing 
section, utilizing the air temperature, wind speed, relative 
humidity, pressure and the cloud cover as forecasted to calculate 

the frost thickness distribution on the wing section as a func- 

tion of time to 6 a.m. By making several computer runs with dif- 
ferent initial atmospheric data at 5 p.m., a graph of the average 
frost thickness versus air temperature, relative humidity, and 
wind speed was derived. Also comparisons were made between frost 
formation on unpainted versus painted wing sections for different 
cloud cover conditions. The frost thickness distribution was 
then converted to lift, drag, and the angle of attack penalties 
at takeoff speed for a general aviation aircraft using the 
empirical formulas from the previous section. 

Applications of the results are as follows; if no aerodyna- 
mic penalties are expected at 6 a.m. as forecasted from 5 p.m.. 
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then no Flight Plan changes are required for the next morning. 

If nocturnal frost with significant aerodynamic penalties are 
expected at 6 a.m., then the pilot can take preventive measures 
such as put the aircraft in the hangar, cover the wings with 
plastic sheets, wait for the sun to dissipate the frost, or can- 
cel the flight. If the preventive measures are not taken, and 
the pilot sees that the wing is covered evenly spanwise with noc- 
turnal frost in the morning, and wishes to continue the flight 
plan, then some modified takeoff procedures are necessary. These 
procedures could be either a reduction in gross weight or an 
increase in the takeoff speed provided, of course, there is suf- 
ficient runway length and the angle of climb is high enough for 
obstacle clearance. An alternative procedure could be the par- 
tial removal of frost from the leading edge to about 10 to 15 
percent of the chord length. This may considerably reduce or 
eliminate aerodynamic penalties. Accidents occur when ordinary 
takeoff procedures are used with a deteriorated wing section. 

This results in an inadequate safety margin above stall speed. 
Lastly, it is to be noted that a takeoff in the morning with noc- 
turnal frost thicknesses uneven across the wing span could also 
be very dangerous. This could produce an asymmetric lift penalty 
across the wind span resulting in a wing drop, or rolloff at 
stall . 
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SECTION 2 


NOCTURNAL FROST FORMATION/DISSIPATION ON A WING SECTION 

In the basic frost formation model by Dietenberger [5] , 
designed for refrigeration type applications, the wall tem- 
perature was held constant as a function of time. On an aluminum 
plate or airfoil exposed to the sky, the plate temperature will 
change with time because nocturnal radiation cooling of the plate 
will produce a plate temperature that is colder than that of the 
surrounding air. Therefore, the plate or wall temperature must 
be calculated as a function of time to determine when it reaches 
the frost-point temperature (Dietenberger [5]) and initiates 
frost. Thus, a model extending the Basic Frost Formation Model 
for calculating the wing section surface temperature due to 
radiation cooling was developed. The input data into the basic 
frost formation model consists of wall temperature, air tem- 
perature, ambient relative humidity, and the heat and mass 
transfer coefficient as a function of time. The model then 
calculates the frost surface temperature, the frost weight, and 
the frost thickness as a function of time. To adapt this model 
to nocturnal cooling, a few changes were necessary. The 
following section discusses the Noctural Frost Formation Model. 

After a frost layer has formed on the aircraft throughout 
the night, the possibility exists that, in the morning, the sun 
may appear prior to the aircraft takeoff. The sun's flux and the 
convective heating may be sufficient to either melt or sublimate 
the frost. In fact, several minutes of sun’s irradiation may be 
enough to remove all evidences of frost. Another source of frost 
dissipation was observed by Langston [1] when a breeze came up 
during the night. Thus, frost dissipation has been incorporated 
into the Noctural Frost Formation Model. 
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2.1 MODELING OF THE V^ING SECTION SKIN TEMPERATURE, T^^^ 

To arrive at a model for calculating the airfoil skin 
temperature, an order of magnitude analysis of heat transfer 
terms was made. This allowed simplification of the heat balance 
equations. There are three heat transfer mechanisms that must be 
considered: radiative cooling of the surface, convective 

transfer of the surface to the air, and thermal conductivity 
within the frost/aluminum/air substances. Two terms that were 
used in this analysis are the Biot number, 

Bi E h WA^I^1„ , (1) 

and the heat content, pCp, of different materials that would show 
a transient temperature property. 

According to the Handbook of Heat Transfer [7] , if Bi^ 0.1 
the temperature distribution within the plate or piece of the 
wing section skin, may be considered uniform with an error less 
than 5%. Consider first the wing section supporting structure. 
Assume it is made of aluminum with a thermal conductivity of 

-i ^,="202 . 4 w/mC. Also assume the conservative estimate of the 
airfoil thickness is approximately W=0.5m. Taking Bi=0.1 as the 
critical value, we can estimate the value of h that would allow 
us to make simplifying assumptions, that is: 

h^_.^ = *^skin = 40.48 W/m^°C (2) 

^ W 

The heat transfer coefficient, h^^, in forced convection is 
typically below 30 Vl/m^°C for a wind speed of 5 m/s. Of course, 
the effective radiative heat transfer coefficient, h^ , introduces 
an additional term for comparison with But h^ is typi- 

cally less than 5 W/m^°C for surface emissivity of 1.0 and at 
freezing temperatures. The equivalent radiation thermal conduc- 
tivity in the air pockets can be estimated by k^ = h^. W, or 

substituting in = 0.2 and in the above values we get 


6 



= 0.5 W/m^°C as the upper limit. The air thermal conduc- 
tivity/ of about 0.025 W/m°C and can be ignored as they are 
much smaller than Thus we need to be concerned about the 

conductive temperature gradient in the wing spar structures along 
thickness of the wing section only in an approximated fashion. 
Also the total mass of the spars is about one-half of the wing's 
mass . 

If the wing aluminum skin thickness is typically about 
W = 0.3 cm, which gives h^j-it=6747 W/m^°C, then the temperature 
gradient across the thickness can be considered to be zero. On 
the other hand, if we take W=chord length of the airfoil of about 
5m, then 4.048 W/m^°C. This is a typical natural convec- 

tion value of h„, which means a strong enough chord-wise distri- 
bution of the temperature gradient will exist on the wing section 
skin. Thus the study of the Biot number shows that the macrosco- 
pic dimensions of the airfoil will show perceptable degree of 
temperature distribution, but in the microscopic dimensions, such 
as that of aluminum slab, the skin temperature can be considered 
to be uniform. 

If any material has a heat content, pCp, and a thermal 
conductivity much lower than other materials, then it can be 
ignored in the analysis. For example, the air pocket within the 
wing section has a heat content, Pa^p = 0.001397 J/m3°c, as com- 
pared to the heat content of the aluminum, *^PA1 ~ 

2.368J/m^°C. Likewise, air has thermal conductivity of 0.25 
W/m°C compared to k^^ = 202.4 W/m°C. The smallness of the air 
thermal conductivity drastically reduces the bulk effect of the 
air pockets. Thus, the air pocket can be neglected as a heat 
sink without introducting serious errors in calculating the wall 
temperatures. On the other hand, the heat content of ice is 
p^Cp^ = 1.88 J/m3°c and a thermal conductivity of about 2.5 
W/m°C. It is noted that the heat content and thermal 
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conductivity of frost would be between that of air and ice, and 
if the frost density is sufficiently high, then the frost cannot 
be ignored as a heat sink without introducing errors in the wall 
temperatures. If, in addition, a fuel tank is located in some 
section of the wing and the tank is filled with kerosene, the 
heat content of the tank would be that of kerosene, PkCpj^ = 

1.52 J/m3°c. Thus kerosene in this case could be an important 
heat sink/source if there is an efficient conduction path to the 
skin. It is noted, however, that the thermal conductivity of 
kerosene, kj^ = .1748 w/m®C, is much lower than that of the alumi- 
num, k^j^ = 202.4 W/m°C. Thus, initially, kerosene would not be 
an important heat sink/source. It is noted if the tank is clo- 
sely attached to, say, the leading edge of the airfoil, and is 
filled with kerosene, then it may be possible to retard frost 
formation overnight by preventing the wall temperature from 
reaching the frost-point temperature. 


Based upon the order of magnitude analysis just provided, 
the formulation of the wall temperature, T„.’, should be con- 
sistent with the following; 


1) The spanwise and chordwise temperature distribution on 
the airfoil skin is a slowly changing function of time 
and position. 

2) The ith aluminum slab area is at uniform temperature 
so that it can be treated as a lump system. 

3) The aluminum slab, the spar supporting structures, and 
the frost layer are significant heat sources/sinks. 

4) The effect of the fuel tank on the airfoil skin tem- 
perature distribution is considered minor due to typical 
locations and conduction path to the surface. Later, it 
is shown how this assumption can be relaxed. 

5) Heat radiation exchanges within the airfoil are negli- 
gible at around freezing temperatures. 
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since a lump system analysis is being pursued, several spe- 
cial cases of heat transfer balances on an aluminum slab are 
considered. The heat balances equations, for example, differ for 
a non-frosted and a frosted slab. If the dew point is above the 
frost point, condensation may occur and then freeze. This 
requires a modification in the heat balance equation. Also the 
attachment or the non-attachment of the supporting structures to 
the aluminum slab will result in different heat balances. Even 
for a frosted slab there are different heat balances between when 
the frost is forming and dissipating. It is our objective to 
treat all possible cases. The following discusses the heat 
balance equation used to obtain the aluminum skin temperature for 
each case. The skin temperature is then used as the wall tem- 
perature in the frost formation model. 


Case A . Clean slab of an airfoil segment not connected to 
supporting structures is shown in Figure 1. Four energy terms 
are involved. They are the enthalpy rate, convective heat, 
radiative heat, and the conductive heat within the skin. Thus we 
obtain the heat flux equation in the difference form in order to 
solve for T^g, as a function of time: 

0a1=pA1«1 dt^ =hHl<’'a-T„.) + »6jltT„'‘(l-E'e) + (%T3‘*+(l-S3)T„tFe 



’■"“Al'^sun+'^Al”! 


“£ - 


'w . 
1 


T T 

Wb - W. 


r 4. I 

^ L.(Lf+L.)/2 Li(L^j+L^)/2 J 


where : 

Jsun “ Ifj(cosfJ cosa^ cos ([)^ + sing sin(j)^) + I^jj(0.42 
+ { [cos (f)^cos 6 - ^JJ2+ [0.58 sin 4 .^] 

+ ( 1- ) Fg ( Ij^sin g+I^g) 

F^ = l/2( 1-sin (j)j^ ) for infinite earth surface area. 


(3) 


(4) 

(5) 
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1/4 


= Tg [1 - (0.5 - 0.006 ✓P^)(1 - K^) J 

= ambient temperature at a height of 1.5 meters 


( 6 ) 



/O .31 

100% 

cirrostratus 



10.63 

100% 

altostratus 



= jo. 85 

100% 

stratus 

(7) 

(0.99 

100% 

nimbostratus 



*0.0 

clear 

night 


Equation 

( 4 ) was 

deve loped from Thermal 

Environmental Engineering 

[ 8] , and 

Equation 

(6) was 

obtained from 

Bagdonas, et.al. [6] . 


For the explanation of each symbol, the reader is referred 
to the List of Symbols. Some observations can be made of the 
heat flux equation. The first is that for a quiet, clear, over- 
night situation, where J =0, and h„ is small, the dominant heat 
flux term for cooling of the aluminum slab to occur must be the 
radiation term. The cooling of the slab is retarded if the air- 
foil surface has a very low long-wave emissivity, e ^ . For 

polished aluminum, the e,, typically has a value of 0.05. For 
unpolished and oxidized aluminum, = 0.2. Since most paints 

have long-wave emissivities around 0.85, a significantly colder 
surface and earlier frost initiation would be expected on a 
painted airfoil. Thus to retard frost formation on an airfoil it 
is important to keep the airfoil cleaned and shined to maintain a 
low value of the long-wave emissivity. In the morning, when the 
sun shines on the airfoil, it is desirable to have maximum 
heating of the aluminum slab to melt or dissipate any frost. 

Therefore, the short-wave emissivity, e. , should have a high 

^-‘-sun 

value to maximize the sun's heat flux to the airfoil. For a 

polished aluminum, e.-, is typically at about 0.2, and for 

^•^sun 

unpolished and oxidized aluminum, = 0.5. If possible, a special 

purpose paint could be utilized to obtain a high value of 

e.. , but a low value of e., . 

"•*-sun "-'•e 
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In the general literature of numerical methods on a system 
of first order differential equations, there are usually two ways 
of solving Equation (3) for One way is to use a predictor- 

corrector integration scheme, and the other way is to use back- 
ward differencing or implicit schemes. From Shampine [9] we find 
that for a high transient or "non-stiff" problem, predictor- 
correctors are appropriate, and if the problem is slow changing 
or "stiff," then implicit schemes should be used. Experi- 
mentation with both type of schemes has shown that Equation (3) 
is definitely a "stiff" problem. By making the following 
approximations : 


dT 


w . 

rr-^ = [T ^ - T ^ 
t '' w . w . 


t- At 


J/At 


( 8 ) 


and ; 



Atj3 ^4^^t 

i 




(9) 


and substituting into Equation (3), we then obtain a linear 
system of equations which was solved by standard methods. It 
turns out that for a "stiff" system a proper choice of a backward 
differencing scheme will be stable no matter what the time incre- 
ment is. Consequently, the time increments are limited by the 
accuracy desired in numerical integration. Thus, Shampine [9] 
and others offer a robust code to solve a system of first-order 
differential equations. However, Equation (3) becomes highly 
nonlinear when condensation, ice, and frost growth, and the 
calculation of heat and mass transfer coefficients are eventually 
included in the model. This requires Equation (3) to be iterated 
at a given time step, a feature not available in the robust codes 
examined. However, we note the present iteration scheme utilized 
is conditionally convergent, that is, only in rare cases do the 
panel temperatures fail to converge when the number of iterations 
is large. 
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Case B 


. An ith slab with spar supporting structures 
attached between the ith slab above and the jth slab below is 
shown in Figure 2. The contribution to Equation (3) for the 
solution of due to the spar support is included as follows. 

In a manner analogous to conduction between adjoining slabs, 
a thermal conduction term was added to the right side of Equation 
(3) for conduction through the spar. This term is; 


%par = )Ch 


1(T - T ) 

spar '^j '^i 


( 10 ) 


A heat capacitance factor was added to the left side of equation 
( 3 ) as ; 


q 

^spar 


capacitance 


f '^spar ^ f ^spar , _ '^i 

^ 2 ^A1 *^pAl dt 


( 11 ) 


We note that the ratio < 1 acts as a normalization heat 

IT X 

flux factor to account for the difference in the attached surface 
area versus the ith slab surface area. The other conduction 
terms in Equation (3) also show a normalization heat flux factor 
with respect to the adjoining slabs. Lastly, a kerosene tank 
could be incorporated in the same way as was done for the spar 
structure. The present model does not include this capability. 


Case C . A water film exists on the ith slab when the wall 
temperature, T^^ , is below the dew-point temperature, but at or 
above the freezing temperature. The water film is assumed to be 
very thin, thus not contributing any heat capacitance effects. 

On the other hand the latent heat flux must be included and is 
added to the right side of Equation (3) at the appropriate time 
in the iteration to converge on the T^^'s as; 

^ liquid latent ~ ^e 
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Airfoil 



Figure 2. 


Non-f rested Slab Area, A^, with a Spar 
Structure Underneath. 
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Case D 


Frost or ice exists on the ith slab when the wall 


temeprature , is below the frost-point temperature and below 

the freezing temperature. If the frost-point temperature 
equation is calculated to be above the freezing temperature, then 
it is likely that ice will exist due to freezing of dew. If the 
frost-point temperature is below the freezing temperature, then 
frost growth is automatically predicted. The ice growth (or 
frost/ice melting) is calculated from the net heat flux at the 
ith slab when the ice/frost surface temperature is set at the 
freezing temperature. If the calculated ice growth flux exceeds 
the convective mass flux, then frost growth is predicted instead. 
The above conclusions concerning frost or ice growth is not 
finalized until Equation (3) has converged with respect to the 
wall temperature distribution, T^^ , at a given time step. Since 
Equation (3) is being solved implicitly and iteratively, a 
feedback mechanism from the ice or frost growth is incorporated. 
Thus Equation (3) is further modified by the following heat flux 
.contribution . The heat capacitance flux of the frost/ice is 
added to the left side of Equation (3) as: 


dT 


^frost/ice capacitance 


w . 

= [Pl c J XI- + Pf.(l-B.)c (13) 


The conductive heat flux due to frost or ice is added to the 
right side of Equation (3) as: 


_ 630 r 4. M _ R ^ 1 r f 

‘^frost/ice conductance “ T^ l^Ii ^ “i i J ) /2 


w . W . 

+ — I 

L.(L^^+L.)/2 J 

The latent heat flux due to frost formation is added 
side of Equation (3) in place of Equation (12) as: 


(14) 

to the right 


^frost latent ^s ‘‘^i^ 


( 15) 
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For the latent heat flux due to ice formation Equation (12) 
is retained, A fusion heat flux is also added to the right side 
of Equation (3) as: 


*^ice fusion ^^s ~ ^e^ 


Finally, changes in the wing section surface emissivities 
occur due to the presence of ice or frost. The blackbody 


emissivity 

of 

the 

surface 

is modeled in Appendix I as; 


^surface 

®it 

^ice 

+ 

( ^Ale 

"ice ^ 

>=Il) J + (1 - Bj) 



X [ 

^ice 

+ 

( "Ale 

"ice^ 

g-10 5(xi. + xs. ) j 

(17) 

where 

^li 

and 

Xs 

. are 
1 

ice and 

frost thickness, respectively. 

in meters. 

and 

Bi 

is 

the frost porosity. The above model of 

the 


surface is based on the assumption that the radiative absorption 
is dominant and the radiation scattering is minimized. Thus when 
Xj^ and Xg^ get large, the surface emissivity approaches the 
emissivity of ice. On the other hand, when Xj^ and Xg^ are very 
small, the surface emissivity approaches that of aluminum 
(providing the wing is unpainted). The porosity weights 
emissivity according to blackbody radiation penetration of the 
frost layer. Likewise, solar absorptivity is formulated as: 


®il-“ice ^ “a 1 °‘ice^ ® ^ij + (1 ) 

X la. + (a,T - a. ) e"^°^''li I 

L ice ' Al ice^ J 


This equation is based on the same assumption about radiation 
absorptions through the frost/ice. 

Given the calculation of the heat and mass transfer coef- 
ficients and of the frost density and thickness, the solution of 
Equation (3), with all the additional terms up to Equation (18), 
can reliably compute the wing section surface temperature 
distribution. With additional modifications, the full wing could 
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be modeled. However, the present model considers only a cross- 
section of the wing. 


2 . 2 MODELING OF THE NOCTURNAL FROST ON ITH SLAB 

In the earlier discussion on the heat content, it was shown 

that the heat content of ice is about the same as that of the 

aluminum. Thus the frost is expected to show some significant 

enthalpy heat rate at high frost densities. This situation 

results in a coupling for the energy equations of the frost and 

of the aluminum slab. The heat flux conducted horizontally from 

the other frosted slabs through the frost layer is considered 

quite minor as compared to the heat flux conducted horizontally 

through the aluminum slab because the thermal conductivity of 

ice, k. =3 W/m°K, is much less than the aluminum thermal con- 
ice 

ductivity, k^^ = 202 V\?/m°K. Thus the heat transfer equation for 
frost formation becomes one-dimensional; 


9Ti g 9T^ 

^fi ‘^pf Tt“ “ 3x t^^i IF" J 


( 19) 


with the frost thermal conductivity from Dietenberger [5] . 
Equation (19) is subjected to the boundary conditions: 


dTi 


x = Xc,. " - Ts.) + hM- L 3 (o>^ - ufe.) + J .3 


“^^i ^idx 

+ - ^e^ + (^e ^ ~ ^e " ^siM^ (20) 



x=0 


3t 


and ; 


( 21 ) 


9T. 

1 


at 

since an analytical 
polynomial function 


I ^=^Si 
solution 
approach 



at 

is not 
toward 


( 22 ) 

possible, a numerical and a 
the problem is pursued. 
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which means the heat flux arriving at the slab through the frost 
structure is the heat flux from the ambient air minus the average 
frost enthalpy heat rate. Physically, this makes sense in that 
the frost layer can be treated almost as a lump system. 

To calculate the frost surface temperature, Tg^, the 
constants a, b, and c from the integration of the heat equation 
for the frost are replaced by their equivalents. The result is; 

Tg . Xg ? 3T„ . 3Tg . 

/ ^ K(T) dT = qg^ Xg^ - pf^ Cpj — g-i- L“JE~^ ^Tt "" 

T«i 

An implicit numerical integration scheme was chosen to 

calculate Tg . for two reasons. The first is that Tg . is 

influenced much more by the changes in and Xg^ than by 

changes in the internal enthalpy rate. Secondly, an order of 

magnitude analysis shows that at frost initiation time the first 

two terms in the above equation are much larger than the third 

term. This means an explicit numerical integration of the term 
3T 

®i is computationally inefficient because the roundoff errors 
3t 3Tg . 

will degrade the values of — rr-^ at the onset of frost formation. 
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Thus the term is approximated by; 




Ts^-^'^J/At, 


(28) 


w 

The parameters , . and T are obtained from Equation (3) 

^ dt w ^ 

with the appropriate additional terms. Then using the same 
numerical techniques described by Dietenberger [5] , the frost 
surface temperature is derived. Finally, the frost density and 
thickness are calculated using the same procedure as in 
Dietenberger [5] . 


The frost surface blackbody emissivity, needed as input 
to Equation (20), is related to the fraction of the frost surface 
in interaction with the radiation as given by; 


^fr 


= (1 - B^) [ 


ice 


+ ( £ 


Ale 


e- ) e 
ice ^ 


-10 


>(XJ ,+x 


Si) 




( 29) 


This equation is similar to Equation (17) which takes into 
account the emissivity of the combined aluminum and frost/ice 
structure. Likewise, the solar absorptivity which replaces 
Ufr in Equation (20) is given by; 

otfr- = (1 - B^) [aiceg + ^ “Alg “ “iceg)® + J. 


Admittedly, the radiation model is idealized because radiation 
heat transfer actually occurs throughout the frost layer rather 
than just at the surface. If, however, we take typical nocturnal 
frost thicknesses of .1 to 10 mm, then according to Equation (29), 
the blackbody emissivity is mostly that of ice, and according to 
Equation (30), solar absorptivity is mostly that of the aluminum 
surface, which seems to agree with observations. 

According to Dietenberger [5] , if the frost surface 
temperature, Tg^, rises to a critical temperature at which the 
water droplets forming on the frost surface do not have time to 
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freeze, but, instead, seep into the frost layer, then all the 
moisture arriving at the frost surface during a time increment 
will be assumed to totally result in the densif ication of the 
frost. Thus the densif ication rate above the critical 
temperature, but at or below the freezing temperature, is 
described as; 


9Pf . 


3t 


■JL = 




- ( Ola - ofe. ) 


(31) 


while the frost thickness remains constant during the next time 
increment. Preliminary results of using Equation (27) for noc- 
turnal frost formation indicates the temperature difference bet- 
ween the wall and the frost surface is typically a fraction of a 
degree. With this kind of small temperature difference, it 
appears the water seepage temperature will never be reached until 
the frost surface temperature is practically at the freezing tem- 
perature. Thus for nocturnal frost formation, the critical tem- 
perature or the water seepage temperature is set at the freezing 
temperature . 


If the frost density finally reaches the ice density via the 
diffusion and the water seepage processes, and if the frost sur- 
face temperature as determined by Equation (27) remains below 
freezing, then the frost is allowed to continue to grow in frost 
thickness, but the frost density is held constant at the ice den- 
sity. That is, we have, for the frost thickness as: 


= 




( - 




dt/ p. 


( 32) 


If by Equation (27) the frost surface temperature goes above the 
freezing temperature, and if the frost density is the same as the 
ice density, then the frost surface will begin melting and liquid 
water will run off. In this case, the frost surface temperature 
will remain at the freezing temperature, the frost density will 
remain a constant (that of ice), and the frost thickness will 


e 
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decrease according to how much liquid water runs off. Thus the 
decrease rate in the frost thickness is calculated by; 


- axg . , 

-5F^Pfl = T- hHi'Ta - Tsj) + h„ J 3 

^ w 

* hid - + (s ’'e + <1 - ''e ’ '^si"] 

Tg . 

- / ^K(T) dT/Xg } (33) 

Twi 

where Tg^ is set equal to 273.16. 

2.3 MODELING NOCTURAL FROST FORMATION ON A FLAT PLATE 

In validating the noctural frost formation model, an airfoil 
geometry was inconvenient for various reasons. Since the testing 
was to be done on a roof of a research building, the blackbody 
radiation from the roof and the nearby shed would have been quite 
difficult to evaluate at various view angles. Secondly, it was 
beyond our existing measurement capabilities to measure local 
wall temperature and local frost weight for different airfoil 
panels. To simplify the problem of experimental validation, yet 
provide useful results, a flat plate lying on a styrofoam box 
inclined at an 15° angle was utilized. Thus the blackbody 
radiation was coming from the sky and the earth, and the thermal- 
couples could be attached underneath the plate. Only a single 
mean frost weight measurement was needed. The flat plate design 
then allowed a simplified analysis of the heat and mass transfer 
coefficients and of the nocturnal frost formation modeling. 
Although the flat plate has only a single blackbody view factor 
and simplified flat plate equations, the main purpose here was to 
verify the physical processes modeled within the frost layer from 
which the frost weight and thickness were calculated. 

When we are not dealing with a complicated geometry such as 
the wing section, but, instead, with just a simple plate at a 
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uniform temperature, there is a very simple method of calculating 
the plate temperature. The heat balance of the plate itself is 
given by; 


Pm Com W = q + ae _ [T't (1 - F ) + [e T 

^Al PaI dt ^ plate « e '• e e 

+ (l-e)T‘+jF - T‘*]+a,^ J 

e 00-' e s ■* plate sun 


(34) 


where ; 


'plate = ®t'ioe ' 'ice J 

It. r , r \ -10x_ I 

“plate " ®l-^ices + l^Alg " Hceg) ® ^ J* 

Combining this equation with Equation (26) we obtain; 


"W r 

- = [q.: 


+ a e 


plate '■ » 


[T*t (1-F ) + (e T‘+ + (1 - e ) T *+ ] F - T *+ I 

Loo 0' oo-'e 


U J. 

+ a 1 I. J "Pc c_. X -rr— 1/ [ P, , c„,. W + c__ X /2 1 ( 35 ) 

plate sun f Pf s 3t ■* '■ A1 PAl *^f Pf s' ■* 


If Equation (35) is combined with Equation (20), then we obtain 

the same results as Equation (3) with the added terms but without' 

thermal conduction terms. Thus an equally good alternative 

dTw 

expression for ^£ -~ to Equation (3) is given by Equation (35) 
for a single plate only. Substitution into Equation (27) from 
Equation (35) and Equation (28) provides complete information 
necessary to calculate the frost surface temperature, and, thus, 
also the frost thickness and density. A predictor-corrector 
scheme was used to integrate Equation (35) to obtain the plate 
temperature as a function of time. To validate the nocturnal 
frost formation model, comparisons were made with experiments on 
a slightly inclined plate exposed to nocturnal frost formation, 
as discussed further in a later section. 
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SECTION 3 


MODELING OF THE HEAT AND MASS TRANSFER COEFFICIENTS 


Given the complicated nature of calculating the nocturnal 
frost formation, it was desired to use the simple yet suf- 
ficiently accurate heat and mass transfer coefficients. The air 
flow regimes over a wing section were examined to see if some 
simple boundary layer technique could be utilized. 


3.1 FORCED CONVECTION 


For a wind flowing in the chordwise and the spanwise direc- 
tion over a wing section, the criterias for boundary layer tran- 
sition and for the boundary layer tripping can be evaluated. 
Jacobs and Sherman [10] show some results for the characteristics 
of an airfoil section as affected by variations in the Reynolds 
number. The various airfoils examined generally retained a lami- 
nar flow behavior up to a chord Reynolds number of about 5x10^. 

If the maximum wind speed for frost formation is taken to be 
5 m/s and the average chord length of a general aviation wing is 
taken as 1.5 m, then the chord Reynolds number is about 5x10^. 
This means the typical wind air flow over the smooth wing in the 
chordwise direction is fully laminar. According to Biguria and 
Wenzel [11] , boundary layer tripping for turbulence induction, 
occurs on a flat plate with a frost layer when the roughness 
Reynold's number: 


Re 


Xs 



> 1200 . 


(36) 


Again, taking the maximum wind speed as 5 m/s, the frost 
thickness would have to be greater than 3.6 mm for boundary layer 
tripping to occur. A typical nocturnal frost thickness is less 
than this, which means the boundary layer tripping is not 
expected on the wing section covered with frost. Thus the 
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boundary layer in the chordwise direction will remain laminar for 
the nocturnal frost formation and a simple procedure for calcu- 
lating the laminar heat tranfer coefficient obtained from Frick 
and McCullough [12] was used. Frick and McCullough validated 
their method on a NACA 65A015 airfoil which is fairly represen- 
tative of many general aviation airfoils. The laminar forced 
convective heat transfer coefficient equation is: 


h = 0.700 k, /<S^ 
X a' L 


and 


( 37) 




\1 ^ ^ “ 1 9* 1 7 fS/c f ^ I 8* 

Re '^U Jo '■U J 


dif J 


(38) 

c “e 

The heat transfer coefficient at the stagnation point is based on 
that of a cylinder and is obtained from the Handbook of Heat 
Transfer [7] as: 


h = 

X 


0.57 c^ 


dU 


e ds 


e 


( 39) 


The potential flow code developed by Bristow [13] was used to 
calculate the potential flow velocities, Ug, around the airfoil. 
This potential flow code allows us to take into account the wind 
blowing from either the front or from the rear direction over the 
airfoil . 


To consider the case of wind blowing from the side of the 
aircraft, or across the wing span, classical results of boundary 
layer analysis for a flat plate were used. The heat transfer 
coefficient is correlated in terms of the non-dimensional Nusselt 
number from the Handbook of Heat Transfer [7] as: 


Nu^ = 0.332 Pr^/^ 


for Re < 3.4 x 10 5 


and : 


Cf/2 


Nu = St Re Pr 

X XX 


(40) 

(41) 


St 


= 0 


.9 {1 + {0. 


2 L. . 0897467 + [^-1 J + In [l + | -1 j J | }( 42 ) 


Cp/2 = 0.185/(log Re ) 
E 1 0 


2* 584 


for Re^ > 3.4 x 10 5. 


(43) 
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The mass transfer coefficient in forced convection is 
obtained from the Chilton-Colburns analogy recommended by 
Dietenberger [5] for laminar and turbulent forced convection over 
a smooth and marginally rough frost surface. The analogy is: 

NUj^ = Shy, Le^/^ (44) 

or: (45) 

where Le is the Lewis number. 

3.2 NATURAL CONVECTION 

When the wind speed is low enough, natural convection will 
dominate the flow on the wing. That is, when the wing section 
surface is cooled below the ambient temperature by nocturnal 
radiation, the air adjoining the wing section is also cooled and 
will start flowing due to the heavier air mass of the cooled air. 
The Nusselt and the Sherwood numbers for a fairly horizontal 
plate (that is, the wing section angle of attack less than 
approximately 5 degrees) are given by: 


Nu^ = 0.27 

(Gr^ 


( 46) 

Sh^ = 0.10 

(Gr^ 


( 47) 


where Gr^ is the Grashof number based on chord length. The 
constant 0.27 was obtained from the Handbook of Heat Transfer 
[7] and the constant 0.10 is a fit to Tajama, et.al. [14] data 
for frost formation on a horizontal flat plate. The natural con- 
vection for an inclined plate (that is, the wing section angle of 
attack greater than approximately 5°) was analyzed by 
Dietenberger [5] and is reproduced in Appendix II for obtaining 
the Nusselt and the Sherwood numbers differently than Equations 
(46) and (47). 
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3.3 MIXED CONVECTION 


When the heat and mass transfer coefficients from both the 
forced and natural convection equations are about equal, the air 
flow is in a mixed convection regime. Churchill [15] recommended 
for the cases of aiding flows where the forced air flow is in the 
same direction as the natural air flow, the mixed convection 
transfer coefficient, or equivalently, the Nusselt and Sherwood 
numbers, are well correlated by: 

*^’^mixed “ ^ ^*^^forced^ ^ ^'^^natural^ ^ ' (48) 

with the similar equation for the Sherwood number. For the 
opposing flows, the mixed convection Nusselt number was corre- 
lated as: 


mixed 


forced 


natural 


(49) 


and in a likewise manner for the Sherwood number. When the 

forced convection flow is perpendicular to the natural convection 

flow, or when the plate is horizontal, the mixed convection 

Nusselt number is taken as the maximum of Nu^ , or 

forced 

NUnatural ® likewise manner for the Sherwood number. 

For comparison of, model results to the nocturnal frost 
experiments on a 15° inclined and 40 cm x 40 cm flat plate, the 
mean Nusselt and Sherwood numbers were derived from Equation 
(40), Appendix II, and Equation (48). 


26 



SECTION 4 


COMPARISON OF THE NOCTURNAL FROST FORMATION PREDICTION 

TO EXPERIMENTS 

To validate the nocturnal frost formation model, fourteen 
sets of experiments, seven of which formed nocturnal frost, were 
conducted by UDRI in the winter of 1981. The experimental 
apparatus consisted of two 40 cm x 40 cm square aluminum plates 
placed on two separate styrofoam boxes. The styrofoam boxes were 
designed to minimize any heat flux to the underneath side of the 
plate. Each plate sat at a 15° inclined angle from the 
horizontal. Fifteen degrees was chosen for the following 
reasons. According to the results by Sparrow, et.al., [16] for 
the forced convection regime, the average Nusselt number for a 
flat plate inclined at an angle of fifteen degrees should be 
equivalent to the average Nusselt number of the same plate in a 
horizontal position. Secondly, the experiments were performed on 
the roof of a building with an observation shed nearby. Thus, it 
was possible to orient the plate away from the building and 
minimize radiation heating from the building itself. No similar 
experiment on a wing section was conducted. 

On one plate three very small thermocouples were attached 
underneath and were used to measure the plate temperature. 
Underneath the other plate was a weight measurement transducer. 
Periodically, the weight of the plate with frost was measured by 
temporarily placing a large styrofoam box over the plate to mini- 
mize the weight errors due to gusty winds. The frost thickness 
was measured on both plates with a penetrometer. The emissivity 
of the plate was found to be 0.2. This value was derived from a 
combination of infrared camera and thermocouple measurements on 
the plate. Along with the observations of plate temperatures, 
frost weight and thickness, the meteorological parameters were 
measured nearby. The wind speed was measured with a sensitive 
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hot-wire probe, the air temperature and relative humidity with a 
hygro-thermograph , and the pressure by a microbarograph. The 
degree of cloud cover and cloud types were estimated by visual 
observation. These meteorological observations are shown in 
Tables I to VII as a function of time. The blackbody sky 
radiation temperature given by Equation (6) is functionally 
related to the air temperature measured 1.5 meters from the 
ground rather than at the building height. Thus a correction was 
made to convert the observed air temperature at the top of the 
building to an equivalent 1.5 meter temperature. Only seven of 
the fourteen experimental data resulted in nocturnal frost for- 
mation. Results of these seven tests are presented in Tables 
I-VII and Figures 3-9. i 

The input data from Tables I-VII was used in the nocturnal 
frost formation model to generate a set of frost prediction 
curves as a function of time. These curves are compared to the 
actual frost measurements in Figures 3-9 . The point symbols are 
the observed values. The open circles represent the measured 
plate temperature, and the solid lines correspond to model pre- 
dictions of plate temperatures. The open triangles correspond to 
the frost weight measurements, and the dashed curves correspond 
to model prediction of the frost weight. Lastly, the open 
squares correspond to the frost thickness measurements, and the 
dot-dashed curves correspond to the model predictions of the 
frost thickness. Frost thickness measurements were not made on 
the first data set. 

Although the thermocouple measurements are accurate within a 
few-tenths, of a degree Celsius, the variability in the wind 
speeds has caused temperature measurement oscillations as high as 
2.6°C and about 1°C on the average. When Figures 3 to 9 are exa- 
mined, the predicted plate temperatures prior to frost initiation 
show agreement with the data within a random error of 1 to 2®C. 
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TABLE I 

Meteorology Data Collected During Frost Formation Experiment #4 
and Converted to Metric Units* 


AIR. 


HOURS 


TIME 
( SEC) 

TEMPERATURE 
(DEG K) 

RELATIVE 

HUMIDITY 

AIR SPEED 
(M/SEC) 

PRESSURE 
( J/M*^3) 

FRACTION OF 
CLOUD COVER 
(CLEAR) 

10:30 

P.M. 

0 

272. 60 

0. 75 

1. 03 

99110 

0, 0 

11: 30 

P,M* 

3600 

273. 16 

0. 78 

1. 03 

99110 

0. 0 

12: 30 

A, M. 

7200 

273. 16 

0. 88 

1. 54 

99110 

0. 0 

1: 30 

A. M. 

10800 

274. 27 

0. 90 

1.03 

99110 

0. 0 

2:30 

A.M. 

14400 

273, 16 

0. 88 

1.03 

99110 

0, 0 

4:00 

A.M. 

19800 

272. 60 

0. 89 

2. 06 

99110 

0. 0 

5:00 

A. M. 

23400 

273. 16 

0. 88 

2. 57 

99110 

0.0 
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TABLE II 


Meteorology Data Collected During Frost Formation Experiment *6 
and Converted to Metric Units* 


"aIh 


HOURS 


TIME 

(SEC) 

TEMPEHATURE 
( DEG K) 

RELATIVE 

HUMIDITY 

AIR SPBEO 
(M/SEC) 

PRESSURE 
( J/M*’*3) 

FRACTION OP 
CLOUD COVER 
(CLEAR) 

11:00 

P.M, 

0 

277. 60 

0. 64 

1. 54 

99750 

0. 0 

2:00 

A.M. 

10800 

275. 38 

0. 70 

0. 69 

99750 

0.0 

3:00 

A.M. 

14400 

275. 38 

0. 70 

0. 56 

99750 

0. 0 

4:30 

A*M. 

19800 

274. 27 

0. 72 

0. 77 

99750 

0. 0 

5:30 

A.M. 

23400 

273. 16 

0. 80 

0. 41 

99750 

0. 0 
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TABLE ril 


Meteorology Data Collected During ?ro3t Formation Experiment #7 
and Converted to Metric Units* 


Tnr 


HOURS 


TIME 
( SEC) 

TEMPERATURE 
(DEG K) 

RELATIVE 

HUMIDITY 

AIR SPEED 
(M/SEC) 

PRESSURE 
( J/M**3) 

FRACTION OF 
CLOUD COVER 
( CIRROSTRATUS 

9:40 

P.M. 

0 

274. 27 

0. 570 

1. 02 

100400 

0, 6 

10:30 

P.M. 

3000 

274. 27 

0. 570 

2. 54 

100370 

0. 6 

12:45 

A.M. 

11100 

273. 16 

0. 680 

1. 27 

100270 

0,0 

2:10 

A.M. 

16200 

273. 16 

0, 670 

2. 54 

100240 

0, 0 

3:00 

A.M* 

19200 

273. 16 

0. 660 

2. 29 

100240 

0. 0 

4:00 

A.M. 

22800 

272. 60 

0. 660 

2. 00 

100240 

0. 0 

5:00 

A.M. 

26400 

272. 60 

0. 690 

1. 70 

100200 

0, 0 

6:00 

A.M, 

30000 

272. 05 

0. 310 

1. 41 

100200 

0. 0 

7:20 


34000 

269, 83 

0. 810 

1. 20 

100170 

0, 0 

7:40 


36000 

269, 83 

0, 810 

1. 20 

100140 

0, 0 
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TABLE IV 


Meteorology Data Collected During Frost Formation Experiment *8 
and Converted to iMetric Units. 


AIR 


HOURS 


TIME 

(SEC) 

TEMPERATURE 
(DBG K) 

RELATIVE 

HUMIDITY 

AXR SPEED 
(M/SEC) 

PRESSURE 
( J/M**3) 

FRACTION OF 
CLOUD COVER 
( STRATUS ) 

9:25 

P.M. 

0 

277.05 

0. 73 

2. 21 

98200 

0,1 

10:55 

P.M. 

4200 

276,49 

0. 80 

2. 21 

98200 

0, 5 

11:45 

P.M. 

3400 

276. 22 

0. 30 

2. 21 

98270 

0.0 

12:45 

A.M. 

12000 

275. 94 

0. 81 

2.01 

98270 

0, 0 

1:45 

A.M. 

15600 

275. 94 

0. 32 

2. 21 

98270 

0. 4 

2:45 

A.M. 

19200 

275. 66 

0. 32 

2. 00 

98270 

0, 7 

3:45 

A.M. 

22800 

275. 38 

0. 82 

1. 75 

98330 

0. 0 

4:45 

A.M. 

26400 

275. 38 

0. 83 

I. 75 

98400 

0. 0 

5:30 

A.M, 

29100 

275. 38 

0, 83 

1. 34 

98400 

0.0 

6:10 

A.M. 

31500 

274. 27 

0. 84 

0, 87 

98470 

0.0 

6:40 

A.M. 

33300 

274. 27 

0. 93 

0. 37 

98500 

0. 0 

7:15 

A.M. 

35400 

273. 99 

0. 93 

0. 41 

98530 

0. 0 

7:40 

A.M, 

36900 

274. 27 

0. 90 

0. 41 

98530 

0,0 
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TABLE V 

Meteorology Data Collected During Frost Formation Experiment 1^9 
and Converted to Metric Qnits. 


HOURS 


TIME 

(SEC) 

TEMPERATURE 
(DEG K) 

RELATIVE 

HUMIDITY 

AIR SPEED 
(M/SEC) 

PRESSURE 

( 

FRACTION OF 
CLOUD COVER 
( CLEAR) 

12: 15 

A.M. 

0 

274, 55 

0. 34 

3, 00 

98930 

0. 0 

1:15 

A.M. 

3600 

274, 27 

0. 36 

3. 09 

98930 

0.0 

2:15 

A.M. 

7200 

273, 16 

0, 88 

2. 50 

99000 

0, 0 

2:45 

A.M. 

9000 

273. 16 

0. 88 

2, 01 

99000 

0, 0 

3: 20 

A.M. 

11100 

273, 16 

0. 38 

2, 00 

99000 

0, 0 

3:50 

A.M. 

12900 

272. 88 

0. 38 

3, 19 

99000 

0, 0 

4: 30 

A.M. 

15300 

272. 33 

0. 88 

2. 38 

99070 

0. 0 

5:15 

A.M. 

18000 

272. 05 

0. 89 

2. 88 

99100 

0. 0 

6:00 

A.M. 

20700 

271. 49 

0, 92 

3. 09 

99170 

0, 0 

6:45 

A.M. 

23400 

271. 49 

0, 92 

2, 01 

99200 

0, 0 

7:15 

A.M, 

25200 

271. 49 

0. 94 

2. 21 

99230 

0, 0 

7: 30 

A.M. 

26100 

271. 49 

0. 92 

1. 75 

99230 

0, 0 
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TABLE VI 


Meteorology Data Collected During Frost Formation Experiment #11 
and Converted to Metric Units. 


AIR 


HOURS 


TIME 

(SEC) 

TEMPERATURE 
(DEG K) 

RELATIVE 

HUMIDITY 

AIR SPEED 
(M/SEC) 

PRESSURE 
( J/M**3) 

FRACTION OF 
CLOUD COVER 
{ STRATUS ) 

10:00 

P.M. 

0 

273. 16 

0. 88 

1.00 

99560 

0. 0 

11:00 

P.M. 

3600 

270. 38 

0.39 

1. 52 

99560 

0.0 

12:00 

A.M. 

7200 

271. 22 

0, 89 

1. 52 

99560 

1. 0 

12:35 

A.M. 

9300 

271. 22 

0. 92 

1.00 

995 94 

0, 0 

1:05 

A.M. 

11100 

270. 94 

0. 93 

0. 81 

99560 

0,0 

1:35 

A.M. 

12900 

270. 10 

0. 93 

1.00 

99560 

0. 4 

2:05 

A. M. 

14700 

269, 33 

0. 95 

1. 50 

99560 

0.0 

2: 35 

A.M, 

16500 

269. 55 

0. 93 

1. 52 

99560 

0. 0 

3:10 

A.M. 

18600 

269. 83 

0, 91 

0, 76 

99560 

0. 0 

4:10 

A.M. 

22200 

269, 27 

0, 39 

0. 51 

99526 

0. 0 

5:10 

A.M. 

25300 

268. 72 

0. 38 

1. 00 

99526 

0, 0 

6:10 

A.M. 

29400 

268. 72 

0. 88 

1.00 

99560 

0. 0 

6:15 

A.M. 

29700 

263. 72 

0. 38 

1. 00 

99560 

0, 0 
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TABLE VII 


Meteorology Data Collected During Froat Formation Experiment #13 
and Converted to Metric Units. 


m 

TIME TEMPERATURE RELATIVE AIR SPEED PRESSURE FRACTION OF 
HOURS (SEC) (DEG K) HUMIDITY (M/SEC) (J/M**3) CLOUD COVER 

(CIRROSTRATUS) 


12: 20 

A.M. 

0 

273. 16 

0, 75 

1. 52 

98781 

0. 0 

1:15 

A.M, 

3300 

272. 05 

0. 75 

1. 52 

98747 

0, 0 

2:00 

A.M, 

6000 

272. 05 

0. 75 

1. 52 

98747 

0.0 

3:00 

A.M, 

9600 

271. 77 

0. 78 

0. 76 

98713 

0, 0 

3:45 

A.M. 

12300 

272. 33 

0. 77 

0. 76 

98713 

0. 0 

4:15 

A.M, 

14100 

272. 60 

0, 78 

0. 71 

98713 

0, 1 

4:50 

A.M. 

16200 

270. 94 

0. 80 

0. 51 

98747 

0. 1 

5:20 

A.M. 

18000 

270. 38 

0, 83 

0. 66 

98747 

0. 1 

6:00 

A.M, 

20400 

270. 66 

0. 82 

0. 46 

98747 

0. 1 

6:30 

A.M. 

22200 

270. 38 

0, 82 

0, 51 

98781 

0. 0 
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Figure 3 . Model Predictions and Data Measurements of Plate Temperature 
and Frost Weight and Thickness as a Function of Time for 
Frost Formation, Experiment #4. 
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Figure 5 . Model Predictions and Data Measurements of Plate Temperature 
and Frost Weight and Thickness as a Function of Time for 
Frost Formation, Experiment #7. 





Figure 7 . Model Predictions and Data Measurements of Plate Temperature 
and Frost Weight and Thickness as a Function of Time for 
Frost Formation, Experiment #9. 








On seven other data sets, not shown, where no frost formed at 
all, similar agreement was found. Once frost starts to form, the 
predicted plate temperatures tend to fall 2.5 to 5.0 degrees 
below the measured values. The more severe cases, (that is, the 
5.0 degrees deviation) generally occurred near sunrise and can 
possibly be associated with solar radiation, which was not 
modeled for these data sets. Nevertheless, even before the onset 
of sunrise, there often exists a 2.5 degrees bias that needs 
explanation . 

On the seven data sets where no frost formed at all, the 
model did not predict any frost formation. In Figures 4, 5, 6, 
and 8, the frost initiation occurred at the correct times as com- 
pared to the data, while in Figures 3, 7, and 9, the frost ini- 
tiation was in error by less than two hours. Thus the frost 
initiation timing is good considering the oscillations in the 
plate temperature caused by the wind speed variability. 

As can be seen from Figures 6 to 9, the accuracy of the 
frost weight and thickness measurements is so good that one can 
identify the shape of the frost weight and thickness curves ver- 
sus time. In all Figures 3-9, except for Figure 8, the predicted 
frost weight tends to increase somewhat too rapidly and the pre- 
dicted frost thickness is somewhat too low. This could be due to 
the predicted plate temperature being too cold since the vapor 
mass flux is in an inverse relationship to the plate temperature. 
The vapor mass flux is also proportional to the heat transfer 
coefficient, but decreasing the heat transfer coefficient would 
produce an even colder predicted wall temperature. If the heat 
transfer coefficient is increased instead so as to increase the 
predicted wall temperature, it may or may not reduce the vapor 
mass flux. There is, however, no other rationale to calculate 
the heat transfer coefficient different than a laminar plate air 
flow assumption. As was demonstrated in an earlier analysis, the 
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boundary layer transition or tripping is unlikely to occur over 
the frost layer to produce a turbulent flow on the plate that 
would give a higher heat transfer coefficient. The possibility 
of laminar boundary layer separation on the blunt styrofoam boxes 
is not feasible because the same wall temperature bias would have 
occurred in the non-f rested case. Thus the presence of frost is 
unlikely to affect the boundary layer over the test plates. 

Allowing the plate wall to see more of the blackbody 
radiation and the frost surface to see less of the radiation in 
the modeling of the surface emissivity will increase the pre- 
dicted wall temperature and decrease the predicted frost weight. 
However, the predicted frost thickness, already too low, would 
get even lower. A higher predicted sky temperature, perhaps due 
to a fog layer, would increase the predicted wall temperature, 
decrease the predicted frost weight, and perhaps increase the 
predicted frost thickness. Indeed, Figure 8 is an exception to 
the rule in providing good predictions of the frost weight and 
thickness which seems to indicate the wall temperature is 
predicted correctly and the surface emissivity model is sound. A 
further point can be made, that if there were systematic errors 
in the heat transfer coefficients or the emissivity model, then 
the same effect on the predicted wall temperature, frost weight 
and thickness as compared to the data should have appeared in all 
nocturnal frost formation experiments. However, the good results 
shown in Figure 8 have provided the exception among the experi- 
ments. An error in the sky temperature allows for exceptions, 
because the earth environment at 1.5 meters height may not corre- 
late well with the temperatures and relative humidity measure- 
ments at building heights in some cases. Somewhat better 
agreement with the data may be obtained if the sky temperature is 
based directly on measurements of air temperature and relative 
humidity at 1.5 m height and the solar radiation taken into 
account. The code for the nocturnal frost formation on the 
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airplane wing as versus that on the plate has incorporated these 
improvements [17] . However, new experiments are needed to verify 
these improvements in the radiation flux to the surface, and to 
uncover any more systematic errors. 

In summary, when the frost weight and thickness prediction 
versus the data are taken together from Figures 6 to 9, we find 
the predictions show the correct qualitative behavior and the 
values of the correct magnitudes. In one instance, in Figure 8, 
the agreement with the measured frost weight and thickness was 
very good. This indicates the blackbody radiation and the heat 
capacitance modeling of the nocturnal frost layer is a good 
approximation to the heat processes involved. As a first known 
model of its kind, the nocturnal frost formation model does quite 
well in spite of the approximations that were made to represent 
the flat plate experiments. 



SECTION 5 


AIRCRAFT TAKEOFF PERFORMANCE ANALYSIS 

This section will evaluate the aerodynamic penalties asso- 
ciated with surface roughness on the aircraft, especially if the 
roughness is caused by nocturnal frost. The aircraft equations 
of motion for a steady climbing flight are then used to analyze 
modified takeoff procedures. Also, an analysis of how a pilot 
not aware of the dangers of frost can or cannot takeoff success- 
fully is discussed in this section. 

Ideally, a potential flow code interacting with a boundary 
layer code that adequately takes into account the presence of 
frost on a multi-element airfoil could be used to evaluate the 
lift and drag coefficients of an aircraft at takeoff speeds. Such 
a program does not presently exist. Thus wind tunnel data on the 
lift and drag penalties due to the upper wing surface roughness 
were examined. Brumby [3] developed some correlations for the 
relative loss in the maximum lift coefficient as a function of 
the average roughness height to chord ratio. Ljungstroem [2] 
measured lift and drag penalties for three configurations of the 
NACA 652A215 wing section as a function of angle of attack, simu- 
lated frost roughness height, Reynold's number, and the roughness 
chordwise distribution. The three configurations consisted of 
the basic wing shape, the wing section with 20° trailing edge 
flaps, and a wing with a high-lift leading edge device. The wing 
with the 20° trailing edge flaps is most representative of a 
general aviation aircraft at takeoff. As far as the actual frost 
roughness height is concerned, Dietenberger [5] has described the 
sand-grain roughness height as equivalent to or up to 66 percent 
greater than the frost thickness. The equivalent sand grain 
roughness generally does not depend on the frost density for a 
typical nocturnal frost layer. 
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The qualitative structure of frost roughness was simulated 
by Ljungstroem using standard aluminum oxide grinding paper. 
Three grades of roughness were used corresponding to roughness 
heights of 0.5, 0.1 and .05 mm. As will be shown in a later 
section, the nocturnal frost formation model predicts a frost 
thickness that tends to be uniform, but covers at times only a 
portion of the wing section. Thus using primarily Ljungstroem 
data, an empirical relationship was established for the drag and 
lift penalties due to the nocturnal frost formation on the wing 
at takeoff speeds. In Appendix III the empirical relationships 
of drag and lift penalties were developed for the NACA 652A215 
wing section with the 20° flaps and a fuselage covered with 
simulated frost. For the same wing section. Figure 10 shows the 
resulting empirical relationship of lift coefficients versus 
angle of attack due to upper surface roughness in comparison to 
L jungstroem ' s wind tunnel data. Figure 11, likewise, shows the 
empirical curves for the drag coefficient versus angle of attack 
due to upper surface roughness in comparison to Ljungstroem wind 
tunnel data. These figures are discussed further in Appendix 
III . 


5.1 TAKEOFF TECHNIQUES WITH ROUGHENED AIRFOIL 

Some researchers, knowing of the consequential takeoff per- 
formance degradation due to the aerodynamic penalties of a 
roughened airfoil, have suggested takeoff modification procedures 
to insure takeoff safety. According to Ljungstroem, "It would 
mean a substantial saving to the airlines if it was sufficient to 
clean only a small part of the wing and control surfaces (that 
is, from the leading edge and 10-20% backwards)." Weeks [4] had 
a significantly different approach: 

"While our knowledge of the detailed flow phenomena which 
are involved is still far from adequate, so that estimates 
such as these can only be tentative ones, it is concluded 
that some combination of take-off weight reduction and 
increase in take-off speed will be necessary in order to 
maintain the normal safety margins in the presence of hoar 
frost deposits." 
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Figure 10. Empirical Curve of Lift Coefficients versus Angle of 
Attack due to Upper Surface Roughness in Comparison t< 


Drag Coefficient 


.06 


.05 


.04 


03 


.02 


.01 L 


0.0 I L_l 

-10 -5 


_l I 1 1 

0 5 10 15 

angle of attack 


Figure 11. Empirical Curves for Drag Coefficients Versus 
Angle of Attack due to Upper Surface Roughness 
in Comparison to Ljungstrem’s Wind Tunnel Data. 
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Finally, Galins and Shirkey's [18] specific recommendation 
for the Boeing 737 aircraft is; 

"When takeoffs are executed during suspected icing con- 
ditions or adverse weather conditions, sound operational 
techniques must be exercised. Wings should be kept clear of 
ice and other forms of contamination and rotation rates 
should not exceed three degrees/second. If deemed 
necessary, improved climb performance (that is, 6% takeoff 
overspeed) may be used for added takeoff speed margins at 
any flap setting." 

The following analysis evaluates each of the suggested modified 
takeoff procedures and, as a result, a new effective takeoff pro- 
cedure is recommended. 

To help evaluate the performance penalties at takeoff, it 
will be assumed the aircraft flies in a steady, linear, climbing 
trajectory after takeoff. It is also assumed that during the 
ground takeoff run, the thrust is much larger than the increased 
drag force due to roughness so that the extra runway length 
covered before takeoff is negligible. An approximated set of 
steady climbing flight equations are: 

Wcos Y = Cj^ {— jsv ( 50) 

P-i 9 

T - Wsiny = (-y JSV^ (51) 

Assume that for any flight at takeoff, the thrust, the air den- 
sity, and the wing area are constants. The lift and drag coef- 
ficients will vary according to the existence of frost on the 
wings and fuselage. The gross weights, the takeoff speed, and 
the trajectory angle are adjustable parameters. Let the 
subscript R stand for rough aircraft surface and the subscript c 
stand for clean aircraft surface. With this situation. Equations 
(50) and (51) are rewritten removing the thrust, density and wing 
area parameters as; 


W cos Y 
c *c 


C V AC 2 

= (1 - 
c 


( 52 ) 


T = 




(53) 


'R 


Equations (52) and (53) relate parameters for steady state 
climbing with a roughened aircraft to those for a steady state 
climb with a clean aircraft. Equation (52), for example, implies 
that if a clean aircraft with weight is climbing at a climb 

angle , then for a roughened aircraft to have a constant climb 

angle of Yd, the Equation (52) relationship must hold. In a 
similar fashion, for an aircraft with a given fixed amount of 
thrust in a steady state climbing conditions, the drag plus 
weight component in the direction opposite to the thrust must 
equal the thrust, for both clean aircraft and roughened aircraft 
configurations. 

Considering the clean aircraft takeoff parameters, W , y , 
Cl^, are known and the roughened aircraft lift and drag 

curves (Cl^^, Cq^^) are available by modifying the clean aircraft 
curves. There remain three unknowns Yj^ and two 

equations. This allows us to assign a value to one of the 
variables and solve for the other two that result in a steady 
climbing flight. 

The analysis was performed on a single engine, general 
aviation type aircraft described by Fink, et al . [19]. The typi- 

cal aircraft parameter values for takeoff are given in Table VIII. 
For a Y(, = 3.5° steady state climb angle, the necessary angle of 
attack and thrust are derived by solving Equations (50) and (51) 
using the lift and drag coeffcients curves versus angle of attack 
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Results ot Takeoff Analysis for a Single Engine, 
General Aviation Aircraft 


































described in Fink, et al . [19] . These values are = 7.47 and 

T = 502.3 lb. These values also correspond to the minimum drag 
to lift ratio of 0.101. 

The next step is to derive a stall margin for a roughened 
aircraft that percentage-wise is the same as for takeoff with a 
clean aircraft. This was accomplished by applying the experimen- 
tal wind tunnel results from Appendix III to the general aviation 
lift curve described by Fink as follows. A clean aircraft will 
take-off at a lift coefficient at a certain percentage less than 
the maximum clean lift coefficient in order to have an adequate 
stall margin. This also implies that the takeoff will occur at 
an angle of attack, , less than the stall angle. Consider that 
the lift coefficients for a roughened aircraft versus angle of 
attack was calculated by taking the clean aircraft (not the 
airfoil) lift coefficient at a clean angle of attack and reducing 
the lift coeffcient by the amount: 

and simultaneously reducing the angle of attack with respect to 
zero lift by the amount: 

Aa = ( Attg/ otg) a, (55) 

which is proportional to We note the relative loss in 

airfoil lift given in Appendix III is a function of 

relative roughness height, (k^/c), percentage coverage of rough- 
ness, P, and the chord Reynolds number, Re^^ , which allows appli- 
cation of the airfoil data to the aircraft and is transparent to 
aspect ratio effects. Thus to obtain the same adequate stall 
margin for a roughened aircraft, or at takeoff a lift coefficient 
at the same certain percentage (that is, 20%) below the roughened 
maximum lift coefficient, will also require a reduction in the 
clean takeoff angle of attack, by the ratio AUg/c^ with respect 
to zero lift. We then get the angle of attack at takeoff for an 
adequate stall margin on a roughened aircraft as: 
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(56) 


otj^ = ( 1 - Aa/a)( - a^) + 

where is angle of attack at the aircraft's (not the airfoil's) 
zero lift (that is, = -4,06). Thus the input to Equations 
(52) and (53) used to evaluate a takeoff procedure modification 
are derived from the above analysis by requiring the same percen- 
tage stall margin as occurred with a normal clean wing takeoff. 

The clean aircraft drag coefficient, for the general 

aviation aircraft as a function of a is obtained from Fink, 

c 

et al. [19]. The roughened aircraft drag coeffcient, in 

Equation (53) is calculated from the clean drag curve as: 

+ aCd^ + = Cd^ + ACq (57) 

where AC^^ is the drag increment on the fuselage due to roughness 
and aCq^ is the drag increment on the wing section due to rough- 
ness. Both AC^j^ and ACp^ are given by empirical relationships 
in Appendix III. Table VIII shows the computed values of bC^/C^, 
Aa/a, ACj^, and at takeoff for the different roughness heights 

0.1, 0.5, and 1.0 mm and the relative coverage of 1.0 and 0.95 on 
the wing section. These values were used as input to Equations 
(52) and (53). 

5.1 TAKEOFF PERFORMANCE EVALUATION 

Galins and Shirkey's suggestion implied the roughened 
aircraft weight to be same as the clean aircraft weight, W = W_. 
Equations (52) and (53) were solved for the climb angle and 
the takeoff speed V„ that sustain a steady climbing flight with 
an adequate stall margin. The results in Table VIII show that 
for roughness in which the loss is aC./C, > 0.3 the climb angle 
becomes negative and the takeoff speed is 1.2 that of the clean 
aircraft. To climb at all for ACj^/Cj^ > 0.3, the stall margin 
will have to be reduced to an unsafe level. Weeks' suggestion 
implies the climb angle remain the same, Yr = and that 
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Equations (52) and (53) be solved for the increased takeoff 
speed, and the relative reduced takeoff weight (W_ - 

Wr,)/W^ that sustain a steady climbing flight with an adequate 
stall margin. The favorable results in Table VIII show that for 
different roughness conditions the takeoff speed, deviates 

from by about ±3%, and the relative reduction in the aircraft 
weight show somewhat the same trends as /C^ . Although this is 
a feasible takeoff procedure modification, it is not the simplest 
approach operationally. Ljungstroem' s suggestion of just par- 
tially cleaning the wings will result in an unsafe takeoff if the 
aircraft is loaded at the maximum clean gross weight, as is 
obvious from Table VIII. The final option, which is the author's 
suggestion, is to maintain the same takeoff speed, V„ = V_ , and 
solve for the climb angle, Yr>/ ^nd the relative aircraft weight 
reduction, (W^ - W^)/W^, from Equations (52) and (53) in order 
to sustain a steady climbing flight with an adequate stall 
margin. The results in Table VIII indicate that the climb angle 
deviates not more than 0.5 degrees from y^, and the required 
relative aircraft weight reduction (W - W„)/W„ is practically 
the same as the relative lift loss, AC^^/Cj^. This means the maxi- 
mum gross weight of the roughened aircraft can be approximated 
as : 

«Rmax " <1 - «cmax <58) 

This takeoff procedure does not have to use explicitly the 
aircraft's equations of motion and is related to partial removal 
of frost through the relative lift loss term, t£. /C . of course, 

Li Li 

while setting V„ = and maintaining constant airspeed during 
climbout, the angle of attack will naturally be at the reduced 
value for a safe margin above stall. 

In this section we present the conceptual conclusions and in 
the last section of this report are the operational recommen- 
dations concerning aircraft takeoff. To summarize the results 
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obtained, the increase in takeoff speed is not always effective 
in compensating for the lift loss and the drag increment due to 
roughened wings. A maximum gross weight reduction, however, will 
always be effective. Thus the most probable cause of takeoff 
accidents with uniform roughened wings, with thrust maintained, 
and with normal takeoff procedures is the loading of an aircraft 
to its maximum allowable gross takeoff weight. Although one of 
the purposes of a safe stall margin is to tolerate small wing 
roughnesses, it is rather obvious the typical nocturnal frost 
thicknesses reported here, and elsewhere can bring the stall 
margin down to an unsafe level (at least for the general aviation 
aircraft examined) when fully loaded. Being at an unsafe stall 
margin can lead to a loss of control in flight when the engine 
thrust is reduced, or when the proper takeoff rotation rate is 
exeeded, or when adequate air speed is not maintained. In this 
type of accident a successful climb is obtained, but a pilot not 
realizing his unsafe stall margin does a "normal" maneuver 
leading to a loss of control of his aircraft. In examining some 
NTSB records, it appears this is the most common type of takeoff 
accident with wing roughness. One example is the Cessna 180 
fatal crash that occurred at the Berlin Municipal Airport in New 
Hampshire in 1977 [20] . The pilot failed to remove the ice and 

snow and took on two passengers. His climb was successful, and a 
normal power reduction after 1500 feet from departure was 
accomplished. Then he lost control and crashed. Another example 
is the crash of a Japan Air Lines McDonnell Douglas DC-8-62F at 
Anchorage, Alaska [21] . The NTSB report reads; "the conditions 
were favorable for the accretion of rime icing (or hoar frost) 
from the time the aircraft approached Anchorage until the crash." 
From the flight recorder data and simulation studies, they found 
the maximum lift coefficient had to be reduced about 15% and 
stall at an angle of attack 2 degrees less than normal for the 
recorded aircraft trajectory. This alone may not have prevented 
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continued flight, but "the aircraft may have been rotated to an 
excessive pitch angle just before it reached takeoff speed, and 
the stall stick shaker did not activate until near impact." 
Finally, "the probable cause of the accident was a stall that 
resulted from the pilot's control inputs aggravated by airframe 
icing while the pilot was under the influence of alcohol." 

In the more severe cases, there is no stall margin 
available, and a pilot would be unaware of it until too late, 
thus leading to an accident. Such appears to be the case in the 
Air Florida Boeing 737 crash in a snowstorm at Washington D.C., 
Jan 13, 1982 [22] . Although wing roughness is due to snow/ice, 
not frost, the conclusions by NTSB are pertinent. There was no 
stall margin because "the aircraft's stall warning stick shaker 
activated almost immediately after liftoff and continued until 
impact." In summary, the aircraft actually had a lower thrust 
than the pilot thought it had and the flight crew apparently had 
only a limited experience with airframe icing and its effects on 
the aircraft. NTSB asserts that "the aircraft could not sustain 
flight because of the combined effects of airframe snow or ice 
contamination which degraded lift and increased drag and the 
lower than normal thrust set by the reference to the erroneous 
EPR indications. Either condition alone should not have pre- 
vented continued flight. Continuation of flight should have been 
possible immediately after stick shaker activation if appropriate 
pitch control had been used and maximum available thrust had been 
added. While the flight crew did add appropriate pitch control, 
it did not add thrust in time to prevent impact." 

A further issue that should be addressed is what the lift 
and drag penalties were as determined from computer simulation 
and flight recorder data, or from a knowledge of wing ice/snow 
thickness as related to the aerodynamic penalties. This 
knowledge of the aerodynamic penalties is needed to compute 
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whether the flight was recoverable or not by application of maxi- 
mum thrust. 

Another issue is that if the pilot knew what his tolerable 
wing roughness for his particular gross weight for a safe stall 
margin was and was able to estimate the roughness height, then he 
would have aborted takeoff. 

To confirm the above conceptual causes of the takeoff acci- 
dents and to test the practicality of the takeoff procedure modi- 
fication, one could suggest three phases of testing. First, if 

necessary, is to re-evaluate the constants in the /Ct and 

^m' ^m 

Aa /a empirical correlations for a specific wing configuration 
s s 

and a specific kind of roughness in a wind tunnel. Perhaps a 
simulator can investigate probabilities of a takeoff accident 
using the algebraic relationships of the type shown in Appendix 
III and using the different takeoff procedures suggested. 

Second, test flight an aircraft with the weight reduction takeoff 
procedure modification using the same kind of roughness and wing 
section. This phase is to test to see if the weight reduction 
takeoff procedure modification is sound for the highly controlled 
wing roughness. Thirdly, take the same aircraft with cleaned 
wings and expose it to the nocturnal frost formation. Measure 
its thickness, percentage coverage, and uniformity of coverage. 
Then, based upon the previous analysis, determine the maximum 
gross weight of the roughened aircraft to assure a safe takeoff 
margin. Finally, takeoff with the aircraft and compare its 
observed performance to that predicted for the observed rough- 
ness . 
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SECTION 6 


PREDICTING OR FORECASTING NOCTURNAL FROST RESULTS 

The nocturnal frost formation model can be run with two 
modes to assess frost growth and aerodynamic penalties. The 
first is the diagnostic mode which uses the observed atmospheric 
profiles such as those in Tables I to VII. Figure 12 is the 
result using the data from Table I at the end of the run. The 
graph shows a wing section temperature and frost density and 
thickness distributions on the NACA 64A215 section at the end 
of the run at 5 A.M. Likewise, Figure 13 corresponds to the data 
input from Table V, and Figure 14 corresponds to the data input 
from Table VI with the predicted temperatures, frost density and 
thickness at the end of the run. Generally speaking, the frost 
on the underside of the wing section has a high frost density 
but a very low frost thickness. The frost on the topside has the 
frost density decrease in the same rate as the mass transfer 
coefficient, but the frost thickness tends to stay at a high 
value. Finally’, the ripple effect in the wing section tem- 
peratures is a reflection of the spar support structures in the 
wing section. 

When Equation (A3. 9) in Appendix III is set equal to zero, 

P, the percentage roughness coverage from the leading edge, equal 

k 

to one, the value of Re ( J at which the aerodynamic penalty be- 
gins to appear is about 0.73. On the other hand, it is recalled 

that on a flat plate the admissible frost thickness roughness 

k 

effect is. Re [ = 100, according to Biguria and Wenzel [11] . 

The pressure gradient on the wing section, in conjunction with 
the roughness effect, may have effectively reduced the admissible 
frost thickness by an order of magnitude. Thus to determine the 
value for P, a conservative approach was taken. That is, 
starting from the leading edge, and moving back, the first frost 
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thickness that has Re ( _> 0.73, will determine P, and then the 

rest of frost on the way back to the trailing edge will be used 
to compute an average frost thickness. 

The averaged frost thickness, P, and a takeoff velocity at 
108 ft/sec is used as an input to Equations (A3. 4) to (A3. 10) in 
Appendix III to determine the drag and lift penalty. Using this 
approach, the percentage maximum lift coefficient lost for the 
Figure 12, frost thickness distribution, is calculated as 16%. 
From the frost thickness distributions as shown in Figure 13, the 
percentage maximum lift loss is 10%; and from the frost thickness 
distributions shown in Figure 14, the percentage maximum lift 
loss is 15%. 


The second mode of the nocturnal frost formation is the 
forecasting of frost at 6 A.M. using the atmospheric data 
measured at 5 P.M. The atmospheric profile is forecasted as 
follows. The pressure, wind speed, and cloud cover are held 
constant throughout the night. The absolute humidity is held 
constant until the air temperature reaches the dew-point tem- 
perature, then the absolute humidity is set equal to the 
saturated absolute humidity at air temperature. The dry air tem- 
perature profile is forecasted as; 


T 

a 


”^o "^o .2ir(t + 11) , 

2 ^ 2 32 J 


( 59) 


2T,-^^ - 0.16853T„ 
rr, _ 5 pm M 

^o 1.83147 


( 60) 


5 pm dj. 

'^M " '^5pm " l-^6am " ^ 2 

taken from Bagdonas , et al. [6]. When the air temperature 
reaches the dew-point, fog or rain is assumed to form. Then the 
air temperature is treated as a constant when it reaches a con- 
densation temperature given by; 
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T 

condense 


minimum [T^ 

^5pm 



12 



(62) 


in degrees Celsius. 

Figure 15 shows the results at the end of a forecasted noc- 
turnal frost formation run. The input data is summarized in the 
figure. Plotted are frost density and thickness distributions on 
the wing section. Also included in the figure are the lift and 
drag penalties. Figures 16 and 17 show the similar types of 
information for other atmospheric conditions at 5 P.M. and fore- 
casted throughout the night until 6 A.M. Figure 17 shows how the 
sun dissipates completely the nocturnal frost, in this case, 
after one and a half hours. By making several more runs for the 
various meteorological profiles, the averaged frost thickness at 
6 A.M. can be plotted as a function of air temperature, relative 
humidity, wind speed, emissivity, cloud cover, and wind direc- 
tion. In Figure 18 are plotted averaged frost thickness versus 
temperature. The wind is from behind the wing in each case, the 
emissivity is 0.2, and the night is clear. The averaged frost 
thickness shows only minor variation with wind, but shows strong 
increases with relative humidity and temperature increase. One 
can draw somewhat similar conclusions for Figure 19 for an 
emissivity of 0.85 (corresponding to a painted surface). V^e note 
the frost thickness can grow as thick as about 1mm, although a 
thickness of 0.5mm is more common. Also shown in Figures 18 and 
19 are the regimes of condensation only and the regime of 
frost/ice growth. The greatest frost thicknesses occur at air 
temperatures of 2.5°C to 5°C, relative humidity of 100%, during 
calm and clear nights. In Figure 18 is shown the open diamond 
symbol corresponding to 50% cloud (cumulus) cover, 70% relative 
humidity, calm winds of 0.5 m/s, and an air temperature of 0°C. 
The frost thickness was about one-third of that of the clear 
night with other conditions the same. In general, wind speeds 
less than 0.5 m/s were clearly in the mixed convection regime, 
but the frost thickness tends to be the same as at 0.5 m/s. 
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Figure 15. A Sample of Nocturnal Frost Forecasting and the Corresponding 
Aerodynamic Penalties, on a NACA 64A215 Wing Section with 
20“ Flaps at 6:00AM. 




o 

O LO 
o 10 
O II 

£ 

II >> O 
4-> LO • 
<1; -r- • O 

i- "O O 

3 *r- II 

•*-> E II 
ro =3 S-. 

S- in "O <D 
O) O) > 

CL O) a; o 
E > 0-0 
OJ •!- CO 
h- -M -a 
ro "O =3 
s- r— c: o 

•I— (U *1— r— 

c q; 3: o 


a. 

Lf> 



o 


VO 


CVJ 


CM 

I 


VO 

I 


o 

I 


(luui) ssau>|o.Lm 


L 

VO 


-L 


J- 


LO 

CD 


O 


cn 


CNJ 


000 

(dd/ 6) /C:^isuaa q.sojj 


O 

CD 


LO 

CVJ 


o 

o 


to 


CVJ 

I 


LO 

I 


10 


I 


o 

CD 


10 

CVI 

I 


4 -> 

to 

O 

o 


“O 

o 


OJ 

> 


td 


O) 


m 0 } 

o -H 


CO 

u 
o 

'D -:j m 
G *w 
(d o CVJ 
5 ^ \ 

fi rH 

O 1 

“H d) iH 

-P Si 


fd Eh 

P 


B 

CD 


u 

-P 


0 • 

M-l 


0) rij 


U} 



-P -H 

P 


CO P 
0 C! 

0 


P D 
Cm W 

0 


iH «S 

ro 


(d 3 

• • 


cor- 


p 0 

D •• 

>1 


+j VD 


0 



0 0) 'O 


;z: rC 

0) 


4J 

-P 



Id 


c 0 

p 


•H +J 

0 


4-> 

O4 


CO (U 

(d 


Id c 

> 


U Q 

P4 


(U 

\ 


p c 'd 


0 0 

(D 


Cm ‘H 

-P 


4J rH 


Lw rd 

0) 


0 

s 

• 

-H 



(U CO 

>1 

c 

iH CO 

pH 

r— 1 

0^-H 

0) 

Cm 

e Q 

4-1 

/d 

0) 

C 

-p 

rH 

d 

ca to 

04 tn 

0 



C P 

0 

m 

< Cm 

u 

0 


<D 

p 

!3 

CJ^ 

•H 

Pm 


Oj,) aunt^Buaduiai aoe^ans uol.:^o^s 6 u|.m 


67 



Averaue Frost Thickness at 6 A.H. 




Average Frost Thickness at 6 A.H. (mm) 



Air Temperature at 5 P.M. (®C) 


Figure 19. Forecasted Averaged Frost Thickness on Top 

of Wing Section at 6:00AM for Surface Emissivity 
of .85 versus Air Temperature, Relative Humidity 
and Wind Speed at 5:00PM. 


SECTION 7 


CONCLUSION 

The research on the nocturnal frost formation on a wing 
section was undertaken to better explain the hazard associated 
with frost on a wing section during a takeoff. This paper has 
provided some pertinent information. A nocturnal frost formation 
model on a wing section has been developed that predicts when the 
nocturnal frost will form and also its thickness and density as a 
function of time. The model shows satisfactory comparison with 
the experimental data on a flat plate. An analysis of the aero- 
dynamic penalties as related to the nocturnal frost formation 
properties was made to determine how much the takeoff performance 
would be degraded by a specific frost layer. With an aircraft 
takeoff assuming equations representing a steady climbing flight, 
it was determined a reduction in the maximum gross weight or a 
partial frost clearance and a reduction in the takeoff angle of 
attack are needed to neutralize drag and lift penalties. This, 
in turn, was related to the nocturnal frost layer. A sensitivity 
study was performed to determine the atmospheric conditions which 
produce the most hazardous frost buildup. 

Under calm conditions, at 100% relative humidity, 5°C air 
temperature, and clear skies, the frost thickness got as high as 
1 mm. On the average, the thicknesses were about 0.5 mm. These 
were the same magnitude of frost thicknesses reported by Weeks 
[4] . For a general aviation aircraft with a takeoff speed of 33 
m/s, this gives a Re(kg/c) ranging from 1204 to 2409. From 
Figure 20, this corresponds to a 33% maximum lift loss, from 
which according to Equation (56) a pilot probably could not make 
a corresponding gross weight reduction. A 20% partial frost 
clearance will only require about a 10% reduction in maximum 
gross weight for a safe takeoff. Using Equations (56) and (58) 
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Figure 20. Relative Lift Loss of the NACA 65A215 Airfoil with 20° 

Flaps (and 25° Slats) versus Surface Roughness Distribution. 


along with the empirical relationship in Appendix III, one can 
evaluate other tradeoff options between partial wing cleaning and 
the maximum gross weight reduction. 

The increase in the takeoff speed was found to be ineffec- 
tive in neutralizing the drag and lift penalties. As a result, 
it was determined that the most probable cause of takeoff acci- 
dents with uniform roughened wings was an aircraft fully loaded 
resulting in unsafe stall margins. The tradeoff performance ana- 
lysis can be applied to most other aircraft, including 
transports, provided the relative lift loss ACl^/Cj^ and the 

relative angle-of-attack reduction Aa/a_ can be properly related 

s s 

to (k /c), P, and Re of a particular aircraft. The drag incre- 
s ^ 

ment due to roughness is only important as to how long it takes 
to reach takeoff speed and to show that an increase in takeoff 
speed was ineffective. Obviously, a wind tunnel and flight 
testing program would need to be conducted to verify the takeoff 
procedure modification before even considering operational appli- 
cation . 


One type of operational application would be to automati- 
cally prevent an inexperienced or an incapacitated pilot from 
taxiing the aircraft if he exceeded his gross weight allowance. 
Also, existing stall warning systems should be modified to take 
into account the occurrence of maximum lift at a reduced angle of 
attack. Of course these applications imply automated measure- 
ments of k and P, and other parameters which is an unlikely near 
term prospect. 

Finally, Figures 18 and 19 show the forecasting of the 
regime of condensation only and the regime of frost/ice growth. 
One could conceivably use these figures to decide on preventive 
measures such as putting the aircraft in the hangar, covering the 
wings during the night with plastic sheets, waiting for the sun 
to dissipate the frost, or cancelling the flight plan. At the 
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At the present, as far as nocturnal frost is concerned, it is 
much more sane and convenient to take the preventive measures 
than to plan a takeoff with frost on the wings and control sur- 
faces , 


73 


APPENDIX I 


CALCULATION OF LONG WAVE LENGTH EMISSIVITY AND SOLAR ABSORPTIVITY 


Expressions for the emissivity and solar absorptivity of 
the wing surface have been developed to account for the presence 
of ice and frost. 


The emissivity for an ice layer is calculated using a 
simple radiative transfer model. The assumptions are that no 
radiation scattering occurs within the ice, and the ice layer is 
uniform in temperature. Radiation is assumed to originate from 
the plate as if no ice layer were present. This radiation is 
attenuated in traversing the ice due to absorption. Thus the net 
radiation from the wing transmitted through the ice layer is: 


ae T 
s s 


(Al.l) 


where ij; is the attenuation coefficient (assumed to be independent 
of wavelength). In addition, each layer of ice of thickness, dz, 
is assumed to radiate as a gray body. This contribution to the 
radiation leaving the upper surface is: 


(1 - R) f^o( ^ dz)T4e-'*'(^-^^ 


( A1 .2) 


where i|; dz = de is the emissivity of the layer dz by use of 
Kirchoff's law and R is the reflectance at the air-ice interface. 
The sum of Equations (Al.l) and (A1.2) is the total radiation flux 
leaving the upper ice surface, and this is equivalent to the 
radiation from a single interface having an emissivity e or: 

q^^j^ = aeT^ = ae^T^e"^^ + ot"^ ( 1 - e"'*^)(l-R) . (A1.3) 

Furthermore, (1 - R) is the emissivity of an infinitely thick ice 
layer through Kirchoff's law and is set equal to 0.94 in the model. 
The value for ip has been chosen as lO^rn"^ [23, 24] . The final 
form of the surface emissivity is found from Equation (A1.3) as: 
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(A1.4) 


e = (Eg - 0.94)exp(-10^x) + 0.94 

The form of the expression for solar absorptivity is iden- 
tical to Equation (A1.4) by use of Kirchoff's law. The solar 
absorptivity of the surface is given by: 

a = ( Og - 0 . 98 ) exp ( -lOx ) + 0.98 . (A1.5) 

For frost, the analysis is more complex because radiation 
scattering can occur, especially at solar wavelengths, and 
because the temperature within the frost layer cannot, in 
general, be considered uniform. The fact that frost is less 
dense than ice implies that less absorption should occur for frost 
than for bulk ice. On the other hand, considerably more scat- 
tering occurs in frost because of the inhomogeneities in refrac- 
tive index caused by its porous structure. This scattering 
results in an increase in path length for radiation propagating 
through the frost layer and a resultant increase in absorption. 
These two effects tend to oppose each other and it is not 
unreasonable to expect both the long wave emissivity and solar 

r 

absorptivity resulting from a thin frost layer to be similar to 
that resulting from an ice layer of the same thickness. 

In the model. Equations (A1.4) and (A1.5) are used to 
determine the emissivity resulting from a thin ice layer on the 
plate surface. In case of frost, these expressions are modified 
to take into account the frost porosity. The modified 
expressions are also used when frost forms over ice. In the 
latter case, the thickness x is the total thickness (ice plus 
frost ) . 
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APPENDIX II 


LAMINAR AND TURBULENT NATURAL CONVECTION ON 
AN INCLINED PLATE 

The Sherwood and Nusselt numbers for Natural Convection are 
reproduced from Dietenberger [5] as modified for an inclined 
plate : 


hwH 

M 


hRH 


Sh„ = = nd+o)^), Nu^= = n/5f (A2.1) , (A2.2) 


H p^D 


where ; 


n = 2/Sc • $ 


Gr 


H 


240 (-|y + Sc $) 


1/4 


(A2.3) 


Gr = ^ 
^ H 2 

V 




J 


T (^1.6453 + 2.6453 ui J 


, (A2.4) 


C = / + 2|| $ + 1 - $ , 


(A2.5) 


= 


1 + (D_ 

o 

1 + 0 ) 


(A2.6) 


The mean Nusselt and Sherwood numbers are: 


NUh 3Nu^ , Shj^ 3^^H* 


(A2.7) , (A2.8) 


NUj^ and Shp^ for turbulent flow are : 


Sh„ = 0.029791 Sc^^^^ 

tl 


Gr 


H 


1+0.49426 <eSc 


2/3 


2/5 


, (A2.9) 




1/3 


(A2.10) 
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t ,sc 2/3 » 

= 1.875 - VO. 765625 + 4.08973 $ ^ .(A2.11) 


The mean Nusselt and Sherwood numbers are: 


Shy 6^'^H ' 


= 6 


(A2.12) , (A2.13) 


77 


APPENDIX III 


EMPIRICAL RELATIONSHIPS FOR THE DRAG, LIFT, 
AND ANGLE OF ATTACK PENALTIES 


To evaluate the minimum drag coefficient of a streamlined 
fuselage covered with roughness at takeoff speed, we note the 
fuselage Reynold's number is usually in the tens of millions for 
a typical general aviation aircraft. It is assumed the fuselage 
is nearly all friction drag and treated by a flat plate model. 
For a smooth flat plate the local friction coefficient from the 
Handbook of Heat Transfer [7] is correlated as: 

-1/7 


Cf /2 = 0.0131 Re, 


for Re > 10 
s 


(A3.1) 


From Young [25] , the correlation of the local friction 
coefficient over a fully rough flat plate is approximated by: 

C£„ = 0.0139 (k for 4 X lO’’ < I"''- <*3.2) 

If P is the fraction of the plate roughened from the trailing 
edge, the relative drag coefficient increase of the plate is 
modeled by integrating Equations (A3*l) and (A3. 2) in the 
expression : 


AC 


D 

'D 


a . 


A 1-P) i ^ r- He 

io ^fs ^ -^(1-P) ^fR 


lo Cfg ds 


- 1 


(A3. 3) 


The result is: 
ACd, 


This equation predicts well the Hoerner's [26] very limited 
data on the minimum fuselage drag due to fully-covered roughness. 
Since Equation (A3. 4) was not fitted to any aircraft data, it has 
some physical significance in that extrapolation is allowed 
beyond Hoerner's data but must remain within the range given by 


[ ^ [ Re 


^ 1/7 _ 


i a 


1 J [ 1 - (1 - p) 


6/7 


j (A3. 4) 
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Equations A3 . 1 and A3 . 2 , Significant nocturnal frost is expected 

to occur only on the upper half of the fuselage. Thus Equation 

A3. 4 was halved to estimate the drag increment, in Table 

VIII for the given values of k , and P was set to one. The 

s 

average frost roughness height, k^, on the fuselage was assumed 
to be same as that on the wing. 

For a plain wing section fully covered with roughness at 
the Reynold's number in the millions, Hoerner's [26] shows a 
correlation of the airfoil minimum drag coefficient versus the 
variable, k^/c, which is the same correlation for a flat plate 
multiplied by a constant. The constant is modeled as a function 
of the airfoil's thickness ratio, and at the 15% thickness ratio 
for the NACA 64 and 65 series airfoil, it shows the minimum drag 
coefficient is nearly all friction drag. Young [25] approximated 
the drag coefficient of a flat plate due to roughness as a one- 

seventh power of the variable ( k /c ) in the range from 4 X 

-7-4 ® 

10 to 5 X 10 from a more complicated expression. Hoerner 

[26] shows data of the minimum , drag coefficients of smooth air- 
foils to be nearly constant for the Reynold's number from about 1 
to 10 million, which is typical for general aviation at takeoff 
speeds. Since nocturnal frost generally will partially cover the 
wing section, it was desired to have a drag increment, 
instead of a drag coefficient so that as the fraction of frost 
coverage goes to zero, the drag increment will go to zero. It 
was found the drag increment approach has another advantage . 

From a cursory examination of Gregory and O'Reilly's [27] , 
Clarius' [28], Ingelman-Ranberg ' s et al . [29], and Abbott and 

Von Doenhoff's [30] data concerning partially or fully roughened 
plain airfoil, it was found the drag difference penalty, 
tends to be a constant as a function of angle of attack until 
near the roughened airfoil stall angle. This assumption of 
constant drag difference would also make the corrections to the 
smooth airfoil drag coefficient due to roughness invariant to 
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aspect ratio effects or other drag components of the aircraft. 
Although the previous observations were for plain airfoils, they 
are assumed to apply also to the NACA 65A215 airfoil with 20° 
flaps partially covered by simulated frost. The empirical fit 
guided by the previous ideas to the Ljungstroem' s data reproduced 
in Figure 11 for the drag coefficients resulted in the equation: 

Cj) = (A3. 5) 

where: = 0.0146 + 1.066 X 10“"^ (a+ 5.9)^, and, (A3. 6) 

ACd^ = [0.1344 (kg/c)^/”^ - 0.02595 ] [l - (l-P)^/^j. (A3. 7) 

It is recommended Equations (A3. 5) to (A3. 7) be used for 
the interpolation of the data. Equation (A3. 7) was used in 
Table VIII to estimate the increment in the wing drag coefficient 
due to the wing surface roughness of the frost layer. The air- 
foil drag coefficient beyond the stall angle is given by [31] as: 

Cp = 2.1 (sina)^*^ (A3. 8) 

In the nocturnal frost formation cases considered in this 
paper, the lift penalty has a severe effect on the climb perfor- 
mance. Brumby [3] correlated the relative maximum lift loss as a 
linear semi-log relationship to the variable kg/c . But when the 
Reynold's number dependences are also considered, a lot of 
scatter was produced in the data versus k /c. When the relative 
maximum lift loss data by Brumby [3] , Ljungstroem [2] , Abbott and 
Von Doenhoff [30], and Gregory and O'Reilly [27] was replotted 

versus the variable Re (k /c ) , on the semi-log scale, the scatter 

c s 

was significantly reduced. Only Ljungstroem data on the NACA 
65A215 airfoil with 20° flaps and with or without 25° slats is 
shown here and is given by Figure 20. The open symbols data are 
without 25° slats and the solid symbols data are with the 25° 
slats. To include the effects of partial coverage of frost, P, a 
relationship similar to Equation (A3. 4) was sought, but having a 


80 


L 


semi-log relationship to the variable. Re (k /c ) resulted in the 

c s 

following empirical relationship as a' fit to the open symbol 
data : 


Ke( k /c ) , / . 

-^ = [.3016 + (.3880 -.3467P)ln JJ [1 -(l-P)-^/^J, (A3. 9) 

Lm 

as given by the solid lines in Figure 20. Equation (A3. 9) should 
be understood for what it is, a reasonable empirical fit to the 
data on the relative maximum lift loss. There were not enough 
solid symbol data to establish any kind of empirical relationship 
for the 25° slats. Thus the dashed lines were simply drawn 
through the solid symbol data. In a slightly, but significantly, 
different manner than Brumby [3] , the relative stall angle loss 
was simply correlated as: 


A a ACr 

-^=Co^ (A3. 10) 

s *-'m 

where the angle of attack is relative to the zero lift point. On 

the data reproduced from Ljungsttoem in Figure 10 for the NACA 

65A215 airfoil with 20° flaps, the value for C was fitted as 

o 

0.725. To generate the curves in Figure 10 in the prediction of 
the lift coefficient as a function of angle of attack due to 
roughness. Equations (54) and (55) in the text, and the clean 
airfoil (instead of aircraft) lift curve was used. That is, take 
any point on the airfoil clean lift curve and reduce the lift 
coefficient value by Equation (54) and simultaneously also reduce 
the angle of attack value by Equation (55). Equations (A3. 9) and 
(A3. 10) were used in Table VIII and Equations (52) and (53) to 
estimate the relative lift loss, and the reduced angle of attack 
for takeoff. 
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